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Entanglement production in coupled chaotic systems is studied with the help of kicked tops. De- 
riving the correct classical map, we have used the reduced Husimi function, the Husimi function of 
the reduced density matrix, to visualize the possible behaviors of a wavepacket. We have studied a 
phase space based measure of the complexity of a state and used random matrix theory (RMT) to 
model the strongly chaotic cases. Extensive numerical studies have been done for the entanglement 
production in coupled kicked tops corresponding to different underlying classical dynamics and dif- 
ferent coupling strengths. An approximate formula, based on RMT, is derived for the entanglement 
production in coupled strongly chaotic systems. This formula, applicable for arbitrary coupling 
fT^ , strengths and also valid for long time, complements and extends significantly recent perturbation 

■ theories for strongly chaotic weakly coupled systems. 
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CN ! I. INTRODUCTION 

(N . 

■ A quantum mechanical system, which consists of at least two interacting subsystems, has an unique property called 
^ \ 'entanglement' yy]. This property is unique in the sense that even if we know the exact state of the system, it is in 

. general not possible to assign any pure state to the subsystems. Entanglement is a nonclassical correlation among the 
' subsystems which exists even between spatially well separated subsystems [3]. This unique property of a quantum 
system has been characterized as a quantum resource for quantum information theory and quantum computation . 
Moreover, quantum entanglement has also been studied extensively from the decoherence point of view. It has been 
\ argued that a quantum system in the presence of an "environment" can loose its coherence and behave more like a 

■ classical system (^]. 

A quantum computer is a collection of many interacting particles. Such a many-particle structure may be prone 
to problems of decoherence and chaos. Decoherence can create some errors in the operation of a quantum computer, 
however, these errors, in principle, can be removed by quantum error correcting codes ^j. On the other hand, 
the problem due to chaos has recently attracted some attention. It has been shown that residual, uncontrolled 

■ interaction between the particles might induce quantum chaos in the quantum computer if the interaction strength 
^ \ crosses certain critical limits and consequently, it may destroy the operational condition of the quantum computer 
O^' ^ ■ Besides quantum chaos can also emerge during the implementation of some quantum algorithms '3| . Obviously, a 

quantum algorithm which simulates a quantum chaotic system is by definition a unitary operation showing quantum 
chaos 0- However, it has been shown that well known algorithms, such as Grover's search algorithm and the quantum 
' Fourier transform algorithm give rise to some unusual combination of quantum signatures of chaos and of integrability 
|6j. The error due to the presence of chaos in a quantum computer can also be corrected by error correcting codes, 
but the presence of chaos enhances the complexity and hence much more error correction is needed Q. Therefore, 
the knowledge of the presence and effects of chaos in a quantum computer is necessary to implement proper error 
correcting codes. Very recently, the behavior of quantum entanglement during the operation of an efficient algorithm 
for quantum chaos have been studied However, here we are interested at a more basic level to study the effect of 
the underlying classical dynamics on entanglement production. 

Recently, several studies have explored this question [lol ITll IT^ ll^ [T3 . ITsl fl^ IT?. 18] . The first one studied the 
entanglement production in an A'^-atom Jaynes-Cummings model 10] and they found that the entanglement rate was 
considerably enhanced if the initial wavepacket was placed in a chaotic region. They also argued that their results 
support an earlier conjecture which predicted that the entanglement production in a chaotic system, coupled to an 
environment, would be more than the regular system |l9j ]. According to that conjecture, the entanglement production 
rate would be higher for a chaotic system coupled to an environment. For the Af-atom Jaynes-Cummings model, each 
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atomic subsystem plays the role of an environment for the other. Later, it has been shown that large entanglement 
production rate is not the hallmark of a nonintegrable system [T]| . Even in the integrable iV-atom Jaynes-Cummings 
model some special initial coherent states exhibit strikingly similar entanglement production as corresponding to the 
chaotic case [l^ . 

In another paper, the entanglement production rate has been related to the classical Lyapunov exponents with 
the help of a coupled kicked tops model They also justified their findings on the basis of the above mentioned 
conjecture [Toj . However, the classical limit of the coupled kicked tops derived in this rather well-quoted work is 
incorrect, in fact it is not even canonical. However they consider very weakly coupled tops and therefore their 
conclusions turn out to be qualitatively valid. In other work, one of us studied the entanglement in coupled standard 
maps and found that entanglement increased with coupling strength, but after a certain magnitude of coupling 
strength corresponding to the emergence of complete chaos, the entanglement saturated pl|. Similar saturation of 
entanglement was also observed for a time evolving state, which was initially unentangled. This saturation value 
depended on the Hilbert space dimension of the participating subsystems, and was less than its maximum possible 
value. It was also pointed out that in analogy with environment induced decoherence, the reduced density matrices 
(RDMs) corresponding to subsystems of fully chaotic systems, are diagonally dominant. 

Later, we derived the saturation value of the entanglement using random matrix theory [l5| . Moreover, we pre- 
sented a universal distribution of the eigenvalues of the RDMs, and demonstrated that this distribution is realized in 
quantized chaotic systems by using the model of coupled kicked tops. Subsequently, an analytical explanation for the 
entanglement production, based on perturbation theory, has been given for two weakly coupled strongly chaotic sys- 
tems [i3 • The authors also found that increase in the strength of chaos does not enhance the entanglement production 
rate for the case of weakly coupled, strongly chaotic, subsystems. In a recent work, entropy production in subsystems 
has been examined as a dynamical criterion for quantum chaos [T^ . It has been observed that the power spectrum of 
the entropy production gets progressively broad banded with a progressive transition from regular to chaotic systems. 
More recently, entanglement production has been investigated in a class of quantum baker's map |l8j |. They also 
found that, in general, the quantum baker's map is a good dynamical system to generate entanglement. 

Besides these studies of entanglement production and decoherence in coupled chaotic systems, extensive studies 
have been done on decoherence of chaotic systems that are coupled to an environment. These studies were mainly 
motivated by the fact that decoherence induces a transition from quantum to classical-like behavior and therefore, 
this decoherent approach can be utilized in a more straightforward way to restore the correspondence between a 
quantum chaotic system and its classical counterpart 19]. Irreversibility is the price of this decoherent model for the 
restoration of quantum-classical correspondence in a quantum system. This irreversibility causes entropy production 
in the system. It has been conjectured, as already mentioned, that this entropy, grows linearly in time with a fixed 
rate determined by the Lyapunov exponents. 

This conjecture has been tested for several model open quantum chaotic systems. It has been shown that the 
entropy production rate, as a function of time, in a quartic double well with harmonic driving coupled to a sea of 
harmonic oscillators has atleast two distinct regimes y^M- For short times this rate is proportional to the system- 
environment coupling strength, and for longer times there is a regime where this rate is determined by the Lyapunov 
exponent. In another work, the entropy production in the baker's map and Harper's map coupled to a diffusive 
environment is studied A regime was found to exist where the entropy production rate is determined by the 

system's dynamical properties like Lyapunov exponents, folding rates, etc., and moreover, in this regime the entropy 
production rate becomes independent of the system-environment coupling strength. Similar results are also reported 
in Refs. [23, HI]. In other work evidence has been presented that the decoherence rate (or entropy production rate) of 
a quantum system coupled to an environment is governed by a quantity which is a measure of both the increasingly 
detailed structure of the quantum distributions (Wigner function) and classical phase space distributions p3 | . 

Very recently, it has been reported that, in open quantum systems, there exists a universal scaling among the 
parameters (effective Planck's constant, measure of the coupling strength between system and environment, classical 
Lyapunov exponents) on which the quantum-classical transition of that system depends P]. In another direction, 
decoherence has been discussed in an open system coupled to a nonlinear environment with finite degrees of freedom 
p6i |. It was found that even though the environment is finite dimensional, the strong nonlinearity of it can destroy 
the quantum coherence. Hence there is a possibility to utilize this finite dimensional chaotic system as a model of 
environment, instead of infinite dimensional heat bath. The above possibility has also been discussed in a recent work 
[27|. Naturally this approach is closely linked to studies like the present one on the coupled kicked tops. 

We have discussed two different approaches in the study of entanglement production and decoherence in chaotic 
systems. First approach was to study the entanglement production and decoherence in coupled chaotic systems by 
performing exact numerical calculation or using some model based on random matrix theory and perturbation theory. 
The second approach was mainly based on approximate master equations. In this paper, following the first approach, 
we have studied entanglement production in coupled kicked tops. We have considered the entanglement production 
for both chaotic and regular cases. Besides considering the effect of different kind of classical dynamics on quantum 
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entanglement, we have also considered the effect of different coupling strengths on entanglement production. We have 
extensively studied a measure of the complexity of the time evolving state, based on the second moment of the Husimi 
function of that state, corresponding to both single and coupled tops. Using RMT, we have explained the behaviors of 
this measure for strongly chaotic cases. We have then derived an analytical formula for the entanglement production 
in coupled strongly chaotic systems using RMT. This formula is applicable for any coupling strength and it also valid 
for sufficiently long time. 

This paper comprises of six sections. In Sec^ we have discussed the quantum and classical properties of coupled 
kicked tops. We have presented the correct classical map of the coupled kicked tops. We have discussed the initial 
states used and have defined the measures of entanglement used here. Finally, we have concluded this section by 
discussing a method to visualize the wavepacket of a coupled system on the phase space of a subsystem. In Sec lIIII 
we have considered a recently proposed method to measure the complexity of a quantum state. Using this method, 
we have defined a measure which quantifies the fraction of the total number of Planck cells occupied by the Husimi 
function of a given state, roughly speaking the amount of "phase space" that is filled by the Husimi. 

The Hilbert space dimension is the number of Planck cell's, each of volume h'^, that fit into the total phase space 
volume. In one-dimension, d = 1, then N — Phase Space Area//i. The above mentioned measure of the complexity of 
quantum states is also approximately equal to the fraction of the Hilbert space occupied by the given state. We have 
observed for the single top that a typical time evolving state can occupy half of the total number of the Planck cells, 
and this happens only for the strongly chaotic cases. Whereas for a highly chaotic top coupled strongly to another 
such top, the above measure, now for the reduced density matrix of each top, has reached a value very close to unity. 
We explain the behavior of this measure, using random matrix theory (RMT), for the strongly chaotic cases. For 
nonchaotic and mixed cases, the time evolving state occupies lesser number of Planck cells and is reflected in smaller 
values of this measure. 

In the next section, Sec lIVI we have presented the numerical results on the entanglement production. In Sec. 
we have derived an approximate formula, based on RMT, to explain the entanglement production in coupled strongly 
chaotic systems. Finally, we summarize in Sec lVIl 



II. COUPLED KICKED TOPS 



A. Quantum Top 

The single kicked top is a system, characterized by an angular momentum vector J = {Jx, Jy, Jz), where these 
components obey the usual angular momentum algebra. We set Planck's constant to unity. The dynamics of the top 
is governed by the Hamiltonian |3l|: 

h 



2j 



n— - 



The first term describes free precession of the top around y— axis with angular frequency p, and the second term is 
due to periodic (5-function kicks. Each such kick results in a torsion about z— axis by an angle proportional to J^, and 
the proportionality factor is a dimensionless constant k/2j. Now, to study the entanglement between two tops, we 
consider the Hamiltonian of the coupled kicked tops which can be written, following Ref. as : 

n{t) = Hlit) + H2{t) + Hi2{t) (2) 
where H,{t) = p,Jy^ + ^ Sit - n), (3) 

Hi2{t) = -J.,J.,^(5(t-n), (4) 

n 

where i = 1,2. Here Hi(t)^s are the Hamiltonians of the individual tops, and Hi2{t) is the coupling between the 
tops using spin-spin interaction term with a coupling strength of e/ j. All these angular momentum operators obey 
standard commutation relations. For the rest of the paper we will only concentrate to the case pi = P2 = 7r/2. 
This special choice of the angular frequencies will simplify both the quantum and classical maps. Since and 
JziS are four mutually commuting operators, the simultaneous eigenvectors of these operators we take as our basis. 
In general, this basis is denoted by mi; j2, "^2) = bij^i) ® |j2,m-2), where Jf\ji,mi) — + l)\ji,mi) and 
Jzi\ji,n^i) — miljiTiTii) . The individual top angular momcntums, ji and j2, could in general be different. 
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The time evolution operator, defined in between two consecutive kicks, corresponding to this coupled Hamiltonian 
is given by, 



where the different terms are given by, 



(5) 



Uf = exp [-i^Jy, 



ut = exp ( -^^Jl I . 



U12 = exp I -i-JziJz2 I , 



(6) 



and as usual z = 1, 2. 



B. Classical Top 

The corresponding classical map of the coupled kicked tops discussed above can be obtained from the quantum 
description with the Heisenberg picture in which the angular momentum operators evolve as: 



J(n + 1) = Ul3{n)UT 



(7) 



Now we have to determine the explicit form of this angular momentum evolution equation for each component of the 
angular momentum. Here we present the time-evolution of Jx^ (see Appendix 



1 



' exp ( -i-Jx2 



■ exp 



J 



Jxi 



{Jzi ~ iJyi) (81 exp ( i-Jx2 



(8) 



The above expression differs from the coupled tops map presented in a previous publication [T^ . Firstly, J^^^ is 
now really a Hermitian operator. Secondly, the terms which arise in the above expression due to the interaction, 
contain J^j^ operator instead of Jy^. We proceed by rescaling the angular momentum operator as {Xi,Yi, Zi) = 
{Jxi 1 Jyt , Jzi )/ ji for i = 1, 2. Components of this rescaled angular momentum vector satisfy the commutation relations, 
[Xi, Yj] = iZi/j, [Yi, Zi] = iXi/j and [Zi, Xi] = iYi/j. Therefore, in j ^ 00 limit, components of this rescaled angular 
momentum vector will commute and become classical c-number variables. Consequently, in this large-j limit, we 
obtain the classical map corresponding to coupled kicked tops as. 



X'2 

Y^ 
Z'2 



Zi cos A12 + 
-Zi sin A12 

Z2 cos A21 + 
-Z2sin A21 
-X2 



Yi sin A12 
f Yi cos A12 

Y2 sin A21 
f Y2 cos A21 



(9a) 
(9b) 
(9c) 
(9d) 
(9e) 
(9f) 



where 



A12 = kXi + eXa; and A21 = kX2 + eXi. 



(10) 



The difference between the map presented above and which was derived in [l^j lies in the form of the angles A12 
and A21. However, these differences are very important. The above map is canonical. It satisfies all Poisson bracket 
relations like {X[^Yl} = Z[,{Yl , Z[} — X[ and {Z[,X[} — Y^ , where i — 1,2 ; and Poisson brackets of any two 
dynamical variables corresponding to different tops are equal to zero. In contrast the classical map presented in [l3l |. 
satisfies the first three Poisson bracket relations, but the Poisson brackets of any two dynamical variables corresponding 
to different tops are nonzero and they are proportional to the coupling strength e, implying that the map is canonical 
only in the uncoupled limit. Moreover, this earlier publication relates the entanglement rate to the sum of the positive 
Lyapunov exponents, which were actually determined using the incorrect classical map. However, they considered 
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FIG. 1: Phase space pictures of the single top, corresponding to different parameter values, are presented, (a) k = 1.0. The 
phase space is mostly regular, (b) k = 2.0. The phase space is still very much regular, but now a thin stochastic layer is visible 
at the separatrix. (c) k — 3.0. The phase space is truly mixed type. Regular elliptic islands are visible inside the chaotic region, 
(d) k — 6.0. The phase space is mostly covered by the chaotic region with few tiny elliptic islands. The solid circle is the point 
at which we will construct the initial wavepacket. 

very weak coupling (e = 10~^) among the tops and therefore the error in the calculation of the Lyapunov exponents 
was very small, these being practically those of the uncoupled tops. Hence we believe that the main conclusions 
presented in that paper are still valid. 

In the limit e — > 0, we will arrive at the map corresponding to the single kicked top, whose Hamiltonian is given in 
Eq. and that map is given by. 



The classical dynamics of the single top have been studied extensively in Ref. j3lL l32l | and is a well studied model 
of quantum chaos. From the above expressions, it is clear that the variables {X, Y, Z) lie on the sphere of radius 
unity, i.e., X"^ + Y"^ + Z'^ = 1. This constraint on the dynamical variables restricted the classical motion to the 
two-dimensional surface of a unit sphere. Following the usual procedure, we can parameterize the dynamical variables 
in terms of the polar angle 9 and the azimuthal angle (j> a.s X = sinScosc/), Y = s'm9 sm cf) and Z — cos 6. In terms of 
this new {9, 0) variables, the above map looks very complicated, and therefore we do not display that map. Moreover, 
during our numerical iterations we use the above three-dimensional form of the map, and after every iteration we get 
back the corresponding {9, (j)) from the relations cos 6* = Z and (f) — tan~^(y/X), where cos 6* and </) are the canonical 
coordinates on the sphere. In Fig. ^ we have presented the phase-space diagrams of the single top for different values 
of the parameter k. For k — \.0 and k = 2.0, the phase-space is mostly occupied by regular orbits. As we further 
increase the value of k, we can see the well known KAM scenario. Finally at k = 6.0, the phase-space is mostly 
covered by the chaotic sea, with very tiny islands. The dark circle, marking the point (6*, (j)) — (0.89, 0.63) in all the 
phase-space diagrams, is representing the point at which we will construct our initial wavepacket. The quantities 
presented in all the figures are dimcnsionless. 
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Z COS kX + Y sin kX 
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X. 
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C. Initial wavepacket 

We use a generalized SU (2) coherent state or the directed angular momentum state as our initial state for the 
individual tops and this state is given in |j, TOj) basis as : 



(j, m,|0o, 0o> = (1 + l7n-^7-''-"' y ( ^ l^^^ ) , (12) 

where 7 = exp(i0o) tan(6'o/2). For the coupled top, we take the initial state as the tensor product of the directed 
angular momentum state corresponding to individual tops. Now on, we will write |j, Wi) as \mi) for notational 
simplification. Explicitly in |m,;) basis this initial product state can be written as [I^ : 

+3 

IV'(O)) = ^ (mi,TO2|V'(0))|mi,m2) 

mi,m2 = -i 

+3 

= ^ {m^\el4>l){m2\el4>l)\mi,m2), (13) 

mi,m2 = -i 

where (mi|6'g, (/)q), i — can be obtained from Eq. (|12ll . 

Now we have the evolution \'4'{'n)) = Urlipin — 1)) = U^\i(j{n — 2)) = .... = U^tlj{0)). Even though, the numerical 
iteration scheme for the above evolutions have already been presented in Ref . |l3l | , here we again present that for the 
sake of completeness. From we have 

{si,S2\'ip{n)) = exp (-1-8182) ^ {si\Ui\mi){82\U2\m2){mi,7n2\il'{n-l)) (14) 

rni,m2 = -j 

where 

(si|C/i|mi) = cxp (-^1^?) (f ) ■ (15) 

rfsiLi (f ) is the Wigner rotation matrix js^l : 



The main problem in calculating the Wigner rotation matrix lies in the calculation of the above sum. Defining that 
sum as Vm, , and starting from V-j — 1 and T^-j+i = 2 si, we can get the other Vmi recursively by using the following 
relation HT^ 



{j - mi + l)V„r,-i - 2siK,i + (j + mi + = 0. 

Besides Wigner rotation matrix can be expressed in terms of Jacobi polynomials and of different Hypergeometric 
functions 35]. However, we have followed the above recursive scheme. 

D. Measures of entanglement 

All the previous studies on the connection between entanglement and chaos, were based on pure states of bipartite 
system, where the von Neumann entropy Sy and the Linear entropy Sr of the reduced density matrices (RDMs) were 
natural measures of entanglement. The definition of these entropies are: 

Sv{n) = -Tri[pi(n)lnpi(n)] = -Tr2[p2("-)lnp2('^)] (17) 
and Snin) = 1 - Tn [pUn)] = 1 - Tr2 [pl{n)] (18) 
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where pi and p2 are the RDMs corresponding to the first and the second top respectively. In the eigenbasis of the 
RDM : 

Sv{n) = -^AJnA, (19) 

i 

SR{n) = l-J^A^ (20) 

i 

where A^'s are the eigenvalues of the RDMs. 

E. Reduced Husimi function 

Since the phase space of the coupled tops is four dimensional (5*^ x S'^), it is not possible to visualize the wavcpacket 
dynamics on such a phase space. Therefore, we use an approximate numerical way to visualize the behavior of the 
time evolving state |'0(n)) in any one of its subspaces. We call this method reduced Husimi function and it is defined 
in the following way. Consider a state {ip) in the angular momentum basis |mi,m2), i.e., 

= ^ amlm2l"^l, m2). (21) 

mi ,m2 

The Husimi function of lip) is \ {zi; Z2|'0)Pj where 

(zi;Z2|V')= X! "'rnim2{zi\mi){z2\Tn2), (22) 
rni ,7712 

and \zi) = \6i,(j)i) are the directed angular momentum states (atomic coherent states). We define reduced Husimi 
function corresponding to first subspace. 



Pih{zi) - / \{zi;z2m'dp{z2), (23) 

where dp,{z2) is the Haar measure : 

dp{z2) = sin92d92d(j)2. (24) 

At: 

Since the phase space of a kicked top is the surface of a sphere of unit radius, the total phase space area is 47r. 
Therefore for the kicked top whose Hilbert space dimension is = 2] + 1, volume of the Planck cell is 47r/(2j + 1). 
Hence the above mentioned Haar measure dp(z) is equal to the number of Planck cells present in the infinitesimal 
area dz — sin 9d9d4>. The integration of dp{z) over whole phase space will give total number of Planck cells N = 2j + l 
present in the whole phase space. One can also write the above expression, Eq. H23|l . as. 



2j + l 



(02>2|^)(V^|^2,02)sinMM02 



(25) 



The above integral is just the partial trace of the density matrix |^/')('(/'| over the second subspace, and hence it gives 
the reduced density matrix (RDM) corresponding to the first subspace. Therefore, 



(26) 



where pi is the RDM of the first subspace. Therefore, the reduced Husimi function is just the Husimi function of 
the RDM. We can write pi = J2iLi K\ei){ei\, where A^'s are the eigenvalues of pi and |ei)'s are the corresponding 
eigenstates. These |ei)'s are also called Schmidt vectors. Therefore, 



N 



Piff(^i,<^i) = 5]A,|(0i,(/.i| 



(27) 
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Thus the reduced Husimi function can also be expressed as the weighted sum of the Husimi functions of the Schmidt 
vectors, where the weight factors are the eigenvalues of the RDM. Identically, we can define reduced Husimi function 
for the second subspace, and is given by, 

N 

P2H(^2,</'2)=5]A,|(02,(/.2|d.)|^ (28) 

where \di)'s are the Schmidt vectors of the second subspace. 



III. SECOND MOMENT OF HUSIMI FUNCTION : A MEASURE OF COMPLEXITY OF A 

QUANTUM STATE 

Reduced Husimi function technique is useful for the visualization of the behavior of the time evolving state on 
the phase space. Moreover, we want a phase space measure of the complexity of any state to relate it with the 
entanglement. There already exists a good measure of that complexity based on the Husimi distribution function, 
Ph = {z\p\z), called 'classical entropy' or Wehrl entropy and that is given by 

SIph) ^ J dn{z)pH\npH (29) 

However, it is difficult to determine the above quantity due to the presence of the logarithmic function. Therefore, 
following a recent proposal we consider inverse of the 'second moment of the Husimi function' W2{ph) as a 
measure complexity of quantum states. This measure is defined as, 

where Nhipn) = j dp{z)PH- (31) 

The quantity W2 represents the effective phase space occupied by the Husimi function of the state p and its unit is 
the Planck's cell volume. We note that a similar kind of quantity, based on the Wigner function, has been introduced 
and studied as a measure of the complexity of quantum states in phase space j30| many years ago. 

We can now define a quantity AiVg^: = W2{ph)/N as the fraction of the total number of Planck cells {N = 2j + I) 
occupied by the state p. Since the total number of Planck cells is equal to the Hilbert space dimension, we can 
define AiVgU also as the rough measure of the fraction of the Hilbert space occupied by the above state. The above 
definitions of AA'gQ- are valid for the single top. For the coupled tops, phase space is 4-dimensional. Here, we can 
define AiVgg for any one of its subspaces. However, the only difference between these two cases is that p is a pure 
state for the single top whereas for the coupled tops, p is a mixed state. Here we have studied the time evolution of 
AA'g-g- for the single top and also for the coupled tops. 



A. Single top 

In the single top case, we have again considered SU{2) coherent state |'0(O)) — |0Oi0o): which we have already 
defined in Ea. (|12|l . as the initial state. We have constructed this state at the point (6*0, ^o) = (0.89, 0.63), and evolved 
it with repeated applications of the single top evolution operator U . The time evolution operator [/, defined between 
two consecutive kicks, is given as 

U = exp [~i\jy) cxp f^-^Yj'^l) ■ (32) 
For the single top case, AA'gg- at time n is 

^^cff = (2j + l)M2[|VXn))] 
where M2[\'^{n))] = [ dfi{z)\{z\^P{n))\'^ (33) 
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FIG. 2: Evolution of AA'gg- is presented for the single top. For the nonchaotic cases { k = 1.0 and k — 2.0 ), denoted 
respectively by solid and dotted line, maximum value of AA'gg- is very less. That means, the time evolving state has very little 
access over the phase space. However, for chaotic cases ( A: = 3.0 and k = 6.0 ), maximum value of AA^'gg- is also not large. For 
the strongly chaotic case { k = 6.0 ), the average value of the maxima is about 0.5. 

and \ip{n)) = U'"'\ip{Q)). In Fig[21 we have shown time evolution of AA^^^g- for different fc-values. For k — 1.0, the 
initial state is inside the elliptic region, and therefore, time evolution of this state is governed by the elliptic orbits 
on which it is initially placed. Since the evolution of this state is in some sense trapped by the elliptic orbits, it has 
little or no access to many parts of the phase space. Consequently, the maximum value of AN^q is very small. After 
reaching its maxima, there are many oscillations in the time evolution of AiVgg due to partial and full revival of the 
time evolving state \tjj{n)). This particular issue of quantum revival of the time evolving state in such mixed systems 
warrant a separate study. Now at A; = 2.0, the initial state is inside a stochastic layer present at the separatrix and 
consequently its dynamics is restricted by and large to be inside that stochastic layer. Naturally, for this case, the 
maxima of AN^q is again small. For k — 3.0, the phase space is of a truly mixed type, with a significant measure 
of chaotic orbits. In this case, the initial state is inside the chaotic region. Therefore, time evolution of this state is 
governed by the chaotic dynamics and this state has access over chaotic region of the phase space. Since the size of 
the chaotic region is large, hence the maxima of AA'gg- is larger (~ 0.35). When k = 6.0, the phase space is mostly 
covered by the chaotic region, with few visible tiny regular islands. The time evolving state has now almost full access 
over the phase space. However, we observed that AA'gg reaches maximum around 0.5 and then fluctuates around 
that value. That means, for this strong chaotic case, the time evolving pure state has access over only half of the 
phase space. This typical behavior of AAgg for strongly chaotic case can be explained by RMT in the following way. 
In the angular momentum basis {|m)}, 

^2[|V'H)] =EE<*IVX'^)>(^(")I^)('I^("))(V'WI™) / dfi{z){z\i){k\z){z\l){m\z). (34) 

i,k l,-m 

After performing the above integral, see Appendix IbI 

M2mn))] = ^^(z|^(n))(V'(n)|fc)(Z|V(n))(V'(n)|m)F(2j;i,fc,Z,m)5,+,,fe+™ (35) 

i,k l.ni 

where F(2rMm) ^ . ) ( ) ( , ) [ ^) i2j ~ ^ - 1)1(2,, + ^ + 1)1 (36) 
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FIG. 3: Distribution of the components of the time evolving state, evolving under strongly chaotic single top dynamics, is 
presented. Top and middle windows are showing that the real and the imaginary part of the components of the time evolving 
state are Gaussian distributed random numbers with zero mean and the variance is l/\/iV, where = 2j + 1 is the Hilbert 
space dimension of the top. In this case j — 80. Bottom window is showing that the distribution of the square of the absolute 
values of the components of the time evolving state are exponentially distributed. This is a typical property of the components 
of a GUE distributed vector. Dotted line representing the GOE (Porter- Thomas) distribution. 



Let us assume, in the angular momentum basis, 

|^(n)) =^c™|m). (37) 

In Figl^l we have presented the distribution of the real and the imaginary part of the coefficients • They are indeed 
Gaussian distributed random numbers. Moreover, in that figure, we have also presented the distribution of |c„ip. 
This figure shows that |cmp are exponentially distributed, which is a typical property of the elements of a Gaussian 
unitary ensemble (GUE) distributed random vector. Therefore, we can assume that the distribution of {c„i} are GUE 
type. For GUE case, RMT average of a quantity identical to M2[|'0('t^)] has been calculated in a recent publication 
|36|. and according to that, 

{M2mn)]) = where N ^ 2j + 1, (38) 

where the angular bracket ( ) represents RMT averaged value. Using the above expression, we have 

iV+ 1 1 / 1 



and for large N limit, 

(ATVgff) ^ 0.5. (40) 
This is the saturation value of AN^q, which was observed in strongly chaotic case k — 6.0. 



B. Coupled tops 



In the last section, we have presented reduced Husimi function technique to visualize the behaviors of the time 
evolving state of the coupled tops on any one of its subspaces. However, to measure the complexity of this state in 
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FIG. 4: Evolution of AA'^gg- corresponding to coupled kicked tops is presented. Solid lines and dotted lines are representing 
the results corresponding to nonchaotic cases { k = 1.0 and k — 2.0 respectively ). Dashed lines are representing the mixed 
case { k = 3.0 ) and dash-dot lines are showing the results corresponding to strongly chaotic case { k = 6.0 ). The top 
window representing the results for the stronger coupling strength ( e = 1Q~^ ), middle window is showing the results for the 
intermediate coupling strength { e — 10"'^ ) and the bottom window is for the weak coupling case { e — 1Q~* ). 



any one of its subspaces, we have to define AA^gg- in a subspace. We have defined AiVg-g- for a given subspace as 



(2j + l)M2(p.ff) 
where ^Vhipin) = J dfi{zi){zi\pi\zi) , 



(41) 
(42) 



and i = 1, 2 representing different subspaces. In FigQJ we have presented the time evolution of the above mentioned 
ANqQ for different dynamics (different k values) and for different couphng strengths e. When couphng strength is 
very weak (e = 10"'*), time evolution of AA'^g- for different dynamics are practically identical to that which we have 
observed in the case of single tops. Therefore, for this coupling strength, effect of the dynamics of one top on the other 
top is very small and two tops are very close to two uncoupled systems. For other coupling strengths, the maxima of 
AA'g-g- has not changed much for the nonchaotic cases (fc — 1.0, and k = 2.0). When e = 10"^, for the chaotic cases 
{k — 3.0, and k — 6.0), AA'gjj first reaches the saturation value which is observed in the case of single tops and then it 
increases approximately linearly with time. However, for the stronger coupling (e = 10~^), it is not possible to divide 
the time evolution of AN^q, for the chaotic cases, into two distinct time regimes. In these cases, AN^q saturates at 
much higher values than the maxima of AN^g observed in single top. For the strongly chaotic case k — 6.0, AA'^g- 
saturates at a value that is very slightly less than unity. This saturation of AiVgg? can be also explained by RMT, 
which we now proceed to do. 

In the angular momentum basis, second moment of the Husimi function of the reduced state, say for the first 
subsystem, at time n is. 



M2{pih) =^^{pl)^kiPl)lm J dp{zi) {zi\i) {k\zi) {zi\l) {m\zi) . 

i.k l.ni 



(43) 



After performing the above integral, we have. 



M2{pih) = ^^{pi)ik{pi)imF{2j; i, k, I, m)Si+i^k+r, 

i,k l,m 



(44) 
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FIG. 5: Distribution of the components of the eigenvectors of the RDM, corresponding to which entanglement production has 
reached the statistical bound. The top and the middle window shows that the real and the imaginary part of the components 
of these eigenvectors of RDM are Gaussian distributed random numbers with zero mean and the variance is 1/VN. Here 
N = 2j + 1 = 161. The bottom window is showing that the distribution of the absolute square of the eigenvectors of the RDM 
are exponentially distributed. Therefore, the eigenvectors of the RDM are GUE distributed. Dotted line representing the GOE 
( Porter- Thomas ) distribution. 

where F{2j; i, fc, I, m) has already been given in Eq. If we write down above expression in the eigenbasis of the 

RDM pi, then we have, 

N 

M2{pih) = ^ XaXp^^{i\(j>a){<l>a\k){^(l)f3){(l)i3\m)F{2j;i,k,l,rn)6i+i^k+rn 

Q(,/3— 1 i,k l,m 

" X! ^« X! X!^*!'^") ('/'a (Z|^a) (0a \m)F{2j; i,k,l, m)6^+i^k+m 

ex. i^k l,m 

+ ^ XgXp y^(z|(/)a)(0a|fc)(/|0/3)((/)/3|TO)f (2j; j, k, I, m)Si+l,k+m 

a, (3 

^ Y.^lQlo.+ E >^o.XpQlp (45) 
where = Y.Y.^\'t>a){<t><.\k){l\(t>o.){(t)c.\m)F{2j;i,k,l,m), (46) 

i,k l^m 

and Qlp = EE<*l<^")<'^"l^)(^l'^/3)('^/3l'^)^(2j;z,fc,Z,m), (47) 

i^k l,m 

where {Xa, \(t>a)} arc the eigenvalues and the eigenvectors of the RDM pi. 

In Fig[5l we have presented the distribution of the real and the imaginary part of the components of the eigenvectors 
{|0q)} of the RDM pi. This figure shows that the real and the imaginary part of {|0a)} are Gaussian distributed 
random numbers. Moreover, FigElalso shows that the distribution of the absolute square of the components of {I^^Iq)} 
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is GUE type. Therefore, from the recent calculation 0|, we can again use RMT average values of and Q\p to 
get RMT average value of M2{pih) as, 

{M2{pih)) 



N - 



N - 



1 



N + 1 



1 



iV+ 1 



(48) 



We know from our earlier work [if 



2N + 1 
7V2 + 2' 



(49) 



Therefore, we have, 



Hence, 



In the large N limit, 



(M2(pm)) = 



1 



Ar + 1 



2N +1 
iV2 + 2 



1 



iV(Af2(pm)) 
(jV+ 1)(A^^ + 2) 
iV(A^2 _^ 2iV + 3) ' 



(A^eff) 



Ar + 1 
iV + 2 



1.0. 



(50) 



(51) 



(52) 



This is the saturation value of AiVgg?, which we have observed in the strongly chaotic (fc = 6.0) and strongly coupled 
(e — lO"'^) case. We emphasize that this is nearly twice that of pure states in a single top. Thus roughly speaking 
the effect of strongly coupling to another chaotic system doubles the phase space access of a state. 



IV. NUMERICAL RESULTS 
A. Classical Phase space 

In FigQ] we have presented the phase space picture of the single kicked top for different parameter values. For 
k — 1.0, as shown in Fig^a), the phase space is mostly covered by regular orbits, without any visible stochastic 
region. Our initial wavepacket, marked by a solid circle at the coordinate (0.89, 0.63), is on the regular elliptic orbits. 
As we further increase the parameter, regular region becomes smaller. Fig^b) is showing the phase space for k = 2.0. 
Still the phase space is mostly covered by the regular region, but now we can observe a thin stochastic layer at the 
separatrix. In this case, the initial wavepacket is on the separatrix. For the change in the parameter value from 
k = 2.0 to fc = 3.0, there is significant change in the phase space. At fc = 3.0, shown in Fig^c), the phase space is 
of a truly mixed type. The size of the chaotic region is now very large with few regular islands. At this parameter 
value, the initial wavepacket is inside the chaotic region. Fig^d) is showing the phase space for k — 6.0. Now the 
phase space is mostly covered by the chaotic region, with very tiny regular islands. Naturally, our initial wavepacket 
is in the chaotic region. 



14 



- j^^l L.^.^[^J J [ .1 I [^ .4. J A^.l |_.^ JL-L L 




2Q0 400 600 BOO 1000 

Time(n) 

FIG. 6: Time evolution of the von Neumann entropy in coupled kicked tops is presented for different coupling strengths and 
for different underlying classical dynamics, (a) e = 10~^. (b) e = 10""^. (c) e = 10~^. Solid line represents k — 1.0, dotted line 
corresponds to fc = 2.0, dashed line is for k — 3.0 and dash-dot line represents k — 6.0. 

B. Time evolution of the quantum entanglement 

In Fig|3 we have presented our results for the entanglement production in coupled kicked tops for the spin j = 80. 
As we go from top to bottom window, coupling strength is decreasing by a factor of ten. Top window corresponds to 
e — 10^^, middle one is showing the results for e = 10^"^ and the bottom window corresponds to the case e = 10^^. 
For each coupling strength, we have studied entanglement production for four different single top parameter values, 
whose corresponding classical phase space picture has already been shown in Fig^ 

1. Coupling e = 10~^ 

Let us first discuss the case of stronger coupling e = 10~^, whose results are presented in Fig|SJi. It shows that there 
exists a saturation of Sy for k — 1.0 and k ~ 2.0, which are much smaller than the saturation value corresponding 
to highly chaotic cases such as when k = 6.0. The saturation value of Sy for k — 6.0 is the statistical bound 
Sv = In(iV) — i ~ 4.57 (where N = 161), which can be understood from random matrix theory p^ . However for 
k = 3.0, corresponding to a mixed classical phase space, Sv is still less than the above mentioned saturation value, 
indicating the influence of the regular regions. 

These two distinct behaviors of the entanglement saturation can be understood from the underlying classical 
dynamics. For k ~ 1.0, the initial unentangled state is the product of the coherent wavepacket placed inside the 
elliptic region [see Fig^a)] of each top. This initially unentangled state will become more and more entangled under 
the repeated application of the coupled top unitary operator Ut- Moreover, if one observes the evolution of the 
reduced Husimi function corresponding to each top, then it can be seen that the initially localized wavepacket starts 
moving along the classical elliptic orbits on which it was initially placed and simultaneously it also spreads along those 
orbits. 
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FIG. 7: Reduced Husimi functions of the time evolving state, evolving under Ut, are presented corresponding to the time at 
which the entanglement production is saturated, (a) k — 1.0. The wavepacket is spread over the elliptic orbits, (b) k = 2.0. 
The wavepacket is spread over the separatrix. It is also showing strong localization at the unstable period-4 orbit, (c) k — 3.0. 
The wavepacket is spread over the whole chaotic region, (d) k — 6.0. At this parameter value, the phase space is mostly covered 
by the chaotic region, see FigH Consequently, the wavepacket is spread over almost whole phase space. 

However, one can observe some initial oscillations in the entanglement production, which is due to the fact that 
the entanglement production is mostly determined by the spreading of the wavepacket along 0-direction. As we know 
cos 9i — linij^aciJzi/ j), therefore the spreading of the wavepacket along ^-direction determines how many eigenstates 
of Jzi, which are also our basis states, are participating to construct the wavepacket. Larger amount of spreading of 
the wavepacket along the 6'-direction causes greater number of basis states to participate in the wavepacket. Moreover, 
coupling between two tops is via interaction between Jj.^ and J^j- Therefore, this interaction term will couple greater 
number of basis states and consequently leads to higher entanglement. 

Initially, the spreading of the wavepacket sometimes may become parallel to the (/)-direction and therefore its 
spreading along ^-direction become less. Consequently, one can observe a dip in the entanglement production. Finally, 
the wavepacket spreads all over the elliptic orbits and the entanglement production reaches its saturating maxima. In 
FigtZ^, we have shown the reduced Husimi function of the wavepacket corresponding to the maxima (saturation) of the 
entanglement production. After reaching its saturation, there are again many dips in the entanglement production. 
These dips are also due to the small spreading of the wavepacket along 6'-direction. However, the localization of the 
wavepacket along ^-direction are now happening due to fractional or full revival of the wavepacket. These revivals are 
actually the single top behaviors which persists even under the interaction with other top. The quantum revivals of 
the wavepacket are interesting phenomena of any quantum system and therefore it requires separate study, especially 
in this rather more complex setting. 

At fc = 2.0, the center of the initial coherent state was inside the separatrix. Therefore, in its time evolution, the 
spreading of the wavepacket was restricted to be inside the separatrix region. Finally, the wavepacket spread over 
the separatrix region, and the entanglement production arrived at its saturation. The corresponding reduced Husimi 
function has been shown in FigEb- Moreover, the reduced Husimi function shows that even though the wavepacket 
has spread over the whole separatrix region, its spread is not uniform. The wavepacket is strongly localized at the 
unstable period-4 orbit. This strong localization of the wavepacket is also a single top behavior which may also 
warrant separate study. 

At A: = 3.0 and k = 6.0, the initial wavepackets were inside the chaotic region. However, the saturation of the 
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FIG. 8: Evolution of the von Neumann entropy, corresponding to the parameter value k = 1.0, are presented for different 
Hilbert space dimensions (TV = 2j + 1). 

entanglement production are different for these two cases. This can be understood as the phase space of the kicked 
top is more mixed type for k = 3.0 than the case k = 6.0. Therefore, the size of the chaotic region is less for k = 3.0 
than its size corresponding to k = 6.0. Consequently, the wavepacket can spread over less of the phase space for 
k = 3.0 than k = 6.0. In Fig[7|:, we have shown the spreading of the wavepacket corresponding to this case. At 
k = 6.0, since the phase space is almost fully chaotic, the wavepacket can spread over almost whole phase space. In 
FigEli, we have shown reduced Husimi function corresponding to this strongly chaotic case. 

As we know, there exists a universal bound on the entanglement for chaotic cases and that bound is given by, for 
the von Neumann entropy, (<S'y)ga,t = ln7iV where 7 = 1/Ve ~ 0.6. Now a natural question is whether there exists 
any such bound on entanglement of the form ln7A^', for the nonchaotic cases like k — 1.0 and k — 2.0. If there exists 
really such an entanglement bound, then what is the N' in terms of TV? We conjecture that TV' is actually the effective 
dimension of the Hilbert space corresponding to those parameter values, i.e., TV' — N^q = ATV^g-TV. Since we know 
the evolution of Sy and of ATV^^g-, we can determine the time evolution of that factor 7 from the relation 

exp{Sv) 1 
^ ~ TV^ ~ TV 

In FiglSland Fig|51 we have shown the evolution of Sy and AN^q corresponding to fc = 1.0 for different Hilbert 
space dimensions. Using the above relation, we determine the evolution of 7 for this fc-value and that is presented in 
FigEl Initially there were some oscillations, later it fluctuates approximately around 7 ~ 0.52 — 0.54 for different 
Hilbert space dimensions. The solid line is showing the average value of 7 at the saturation region. Figs 1 1 II and 1121 are 
similarly showing the evolution of Sy and of ATVgg at A; = 2.0 corresponding to different Hilbert space dimensions. 
In Fig ll3l we have shown the evolution of 7 for this case. This figure is showing that at the saturation 7 ~ 0.40 — 0.43 
for different Hilbert space dimensions. At the saturation, the factor 7 is different for k — 1.0 and k — 2.0. This is 
essentially due to the fact that at fc = 1.0 and k = 2.0, two different kind of dynamics are responsible for the spreading 
of the wavepacket on phase space. At fc = 1.0, the wavepacket has spread over the regular elliptic orbits, whereas 
at A; = 2.0 the wavepacket has spread over a thin stochastic layer present at the separatrix. Even though we may 
not expect any universality in the case of integrable or near-integrable cases, we have found that for a given coupling 
strength and for a given classical dynamical behavior, the factor 7 is more or less same for different Hilbert space 
dimensions. 




AN, 



pff 



(53) 
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FIG. 9: Evolution of AA^'gg-, corresponding to fc = 1.0, are presented for different Hilbert space dimensions. 
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FIG. 10: Evolution of the factor 7 are presented for different Hilbert space dimensions. This factor has been calculated 
numerically using Ea.H58|l. Here k = 1.0. 

2. Coupling e = 10"^ 

Entanglement production corresponding to this coupling strength has been presented in Figdb). For the nonchaotic 
cases {k = 1.0 and k = 2.0), the saturation value of the entanglement production is less than the entanglement 
saturation observed in the stronger coupling case (e = 10"^). For weaker coupling, the influence of one subsystem 
on the other subsystem becomes less, and the individual subsystems behave more like isolated quantum systems. 
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FIG. 11: Evolution of the von Neumann entropy, corresponding to the parameter value k = 2.0, are presented for different 
Hilbert space dimensions. 
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FIG. 12: Evolution of AA'gjj, corresponding to fc = 2.0, are presented for different Hilbert space dimensions. 



Consequently, pure quantum effects play dominant role in the evolution of the wavepacket. In Fig^J we have shown 
reduced Husimi function for k = 1.0 and k = 2.0 at the time n = 384 when the entanglement production saturated. 
For k = 1.0, the reduced Husimi function is showing that the wavepacket has spread over the elliptic orbits, but not 
uniformly. 

Now for k — 2.0, at the entanglement saturation, the wavepacket has spread as usual over the whole separatrix 
region. Moreover, it also shows localization at the same unstable period-4 orbit. However, the difference is that the 
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FIG. 13: Evolution of the factor 7 are presented for different Hilbert space dimensions. This factor has been calculated 
numerically using Eg. 15311 . Here k = 2.0. 



3 

2.5 

■ J 

1 



I 




-3 -: 



"0 — r 



5 



1. 

1 

r . s 

3 




-3 



1 2 J 



FIG. 14: Reduced Husimi functions of the time evolving wavepacket are presented corresponding to the time n = 384 at which 
the entanglement production gets saturated, (a) k = 1.0. The wavepacket is spread over the elliptic orbits, but the spreading is 
not uniform, (b) k = 2.0. The wavepacket is spread over the separatrix and shows strong localization on the unstable period-4 
orbit. Here e = 10^^. 



wavepacket is now more localized at a particular periodic point of that period-4 orbit which was very close to the 
initial wavepacket. As we have seen in Fig. within our observational time {n = 1000), AA'^g- has not reached any 
saturation value for the mixed and as well as for the chaotic cases. Moreover, for the strong chaos case, k = 6.0, the 
ANqQ was well short of unity even after the observational time and consequently the wavepacket has not got access 
over whole Hilbert space within this time of observation. Therefore, the entanglement production is well short of the 
known statistical bound In(A^) — ^. 



3. Coupling e = 10 



The entanglement production for this very weak coupling regime has been presented in FiglJfc). The entanglement 
production for this weak coupling has recently been explained by perturbation theory [T^ . However, the formula for 
the entanglement production presented in that work is not valid for arbitrarily long times. In the next section we have 
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presented an approximate formula for the entanglement production in coupled strongly chaotic systems. This formula 
explains the entanglement production for the case k = 6.0. Here we have also observed that entanglement production 
is much larger for the nonchaotic cases than the chaotic cases. Rather, we can say that, for weakly coupled cases, the 
presence of chaos actually suppresses entanglement production. 



V. ENTANGLEMENT PRODUCTION IN COUPLED STRONGLY CHAOTIC SYSTEM 



Due to the relatively simple form of Sr, the linear entropy, it is easier to derive an approximate formula for its time 
evolution. Here we present an analytical formalism for the time evolution of Sr in coupled strongly chaotic systems. 

Let us assume, the initial state is a product state, given as |V'(0)) = 101(0)) \(j)2{0)), where \(j)i{0))^s are the states 
corresponding to individual subsystems. In general, the time evolution operator of a coupled system is of the form 
U = UJJq = U^{Ui (Si U2), where is the coupling time evolution operator and J7j's are the time evolution operators 
of the individual subsystems. Furthermore, we have assumed 



{/e = exp(-iei?i2) 



(54) 



where H12 = h'^^^ (g) h^'^\ and the h^^^ are Hermitian local operators. For simplicity, we derive our formalism in the 
eigenbasis of /i^*^'s, i.e., h^^^\ea^) = edea^), where {ea\ \sa^)} are the eigenvalues and the corresponding eigenvectors 
of 

The one step operation of U on |V'(0)) will give, 



(55) 



where IV'(I)) is the time evolving state of the full coupled system at time n = I and |'0o(l)) is the same for the 
uncoupled system. From the above expression, one can get the RDM corresponding to one subsystem by tracing over 
the other subsystem. The RDM corresponding to first subsystem is given by. 



[pi(l)W = (e«h(l)|e«) = ^(eW,ef|V'(l))(^(l)|e«,ef 



e«,ef|^o(l))(v.o(l)|e«,ef\. 



(56) 



Here we now assume that, |'!/'o(l)) is a random vector. Consequently we can further assume that the components of 
|'i/'o(l)) are uncorrelated to the exponential term coming due to the coupling. Hence we have, 



^ ^ E (e«, |V;o(l)) (^o(l)|e(^\ e(2)) Eexp [-ie (e« - e«) 

7 5 

= ^[pio(l)]a/3 E«^P ' - 4'^] ' 



=,(2) 



(57) 



where N is the Hilbert space dimension of the first subsystem and pio is the density matrix corresponding to the 
uncoupled top. If we proceed one more time step, then at the time n = 2 we have. 



[pi(2)]a/3 ^ ^|p(e)l'[pio(2)]„0E°^P 

7 

where p{e) = -^J2^'^p{~^^^^'^^^?)- 



) _ e(i) 



If we use the above assumptions upto any arbitrary time n, we obtain 

1 



[Pi{n)U = ^b(e)|^("-^)[pio(n)]„;3Ee^P ['^^ (^i'^ ' ^ 



(58) 



(59) 



From the above expression, it is straightforward to calculate Linear entropy and that is given as, 

^ 1 - ]^ bwr^"-^^ E E 

a,p 7,(5 



-^e(eW-eW) (e(^) - ) 



(60) 
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FIG. 15: Evolution of the Linear entropy for the coupled strongly chaotic system is presented. The dotted line is the numerical 
results of the coupled kicked tops system. We choose k = 6.0 for the first top and A: = 6.1 for the second top. The solid line is 
the theoretical estimation, given by Ea. lf?)^ . 



This is a general result, applicable to any coupled strongly chaotic systems of the form Ue{Ui (g) U2)- Moreover, this 
result is valid for long time. 

For the coupled kicked tops H12 — Jzi ® Jz^l ]■ Therefore, for this particular system, the above formula would 
become, 



exp 



-i-{mi - ni){m2 - ^2) 
J 



where 



1 ( t \ 

P(^) ~ E exp f — i-mim2 j and = 2 j + 1 



mi,m2 = -J 



(61) 



In large j-limit, we can substitute above sums by approximate integrals and then performing those integrals we get 
(for details, see Appendix 0, 



5«(n)^l-p(e)4("-i) 



2 f Si(2iVe) 
iV 1 e 



1 



{1 - cos(2A^e) + Ci(2A^e) - ln(27Ve) - 7} 



where p(e) ~ — 



(62) 



The functions Si and Ci are the standard Sine-integral and Cos-integral function, while 7 = 0.577216 ... is the Euler 
constant. In the above derivation wc have not assumed, unlike the perturbation theory [l^ . any particular order 
of magnitude of the coupling strength e. Therefore, as we demonstrate below the above formula is applicable for 
non-perturbative coupling strengths as well. 

In Fig ll5l we have shown the numerical result of the Linear entropy (Sr) production in the coupled tops where 
the individual tops are strongly chaotic. Here we have considered many initial coherent states at different parts of 
the phase space and presented the Linear entropy production, averaged over all these initial states, with time. In 
all our previous calculations we only considered the entanglement production on coupling identical tops, therefore, 
permutation symmetry was present. As in the above derivation, we have not assumed any special symmetry property, 
we break permutation symmetry by taking slightly non- identical tops with k — 6.0 for the first top and fc = 6.1 for 
the second top. 
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Fig^] demonstrates that our theoretical estimation, denoted by the sohd curve, is not only valid for weak coupling 
case like e = 10"'' but it also valid for sufficiently strong coupling cases like e — 10~^. Moreover, this formula is 
applicable for very long times. If we consider weak coupling approximation, i.e., je ^ 1, then the above formula will 
become approximately, 

+ (63) 

Therefore, at this weak coupling approximation, the entanglement production rate is 2e^j^/9, which has been calcu- 
lated in a recent publication jj^] by very different means. 

VI. SUMMARY 

In this paper, our major goal was to study entanglement production in coupled kicked tops. Single kicked top is a 
well studied model of both classical and quantum chaotic system. The classical map corresponding to coupled kicked 
tops was presented in a previous publication, but was unfortunately incorrect. Hence, we have presented the correct 
classical map corresponding to the coupled kicked tops which is canonical. In the quantum case, we have studied the 
reduced Husimi function to visualize the behavior of the wavepacket of a coupled system in any one of it subspaces. 
We have also studied a phase space based measure of complexity of the time evolving state ( denoted by AiVgjj ), 
which quantify the fraction of the total number of the Planck cells occupied by the Husimi function of a given state. 
As we have already mentioned that, for kicked top, this quantity is also approximately equal to the fraction of the 
Hilbert space occupied by a given state. We have studied this quantity for both single and coupled tops. It has been 
observed that, for the single top, the time evolving state can occupy maximum, in average, half oi the total number 
of the Planck cell, i.e., AA'g-g: — 0.5, and this happened for the strongly chaotic cases. 

For nonchaotic and mixed cases, the time evolving state occupies even less number of Planck cells and it is reflected 
in smaller values of AN^q. Following a recent result, using RMT, we have explained the fact that AN^g = 0.5 for 
the time evolving state corresponding to strongly chaotic single top. However, when a strongly chaotic top is strongly 
coupled to another such top, AN^q corresponding to any subsystem reaches very close to unity. We have again 
explained this by means of RMT calculations. 

Then we studied entanglement production in coupled kicked tops for different underlying classical dynamics of 
the individual top and also for different coupling strengths. We find, in general, entanglement production is higher 
for stronger chaotic cases. Moreover, coupling strength between two tops is also an important parameter for the 
entanglement production. For example, when the coupling strength between two tops is very weak, we find that 
entanglement production is higher for sufhciently long time corresponding to nonchaotic cases. Finally, we have 
derived an approximate formula, based on the ideas of RMT, for the entanglement production in coupled strongly 
chaotic system. This formula is applicable, unlike perturbation theory, to large coupling strengths and is valid for 
sufficiently long times. 

APPENDIX A: DERIVATION OF EQ.© 

Let us define ladder operators, 

Ji± = ± Jj/i ; Ji+ = Ji- 
Ji+|mi) = C„Jtoi + 1) and Ji_ |mi) = D^Jmi - 1) (Al) 

where Cmi and -Dmi are known functions of j and mi and |mi) arc the standard angular momentum basis states. 
We can write Jx^ = + Ji-)/2. Therefore, 

4, lul{Ji+ ® hWr + ^4(^1- ® I2)Ut, (A2) 

where the terms present at the right hand side are the Hermitian conjugate of each other. Therefore, it is sufficient 
to determine only one term. Here we will calculate the first term explicitly. We have 

4(Ji+ ® h) = (C/i ® U2)^uf^{Ji+ ® hM^iUi ® U2). (A3) 
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In |mi,m2) basis, Ui2{Ji+ <S) -^2)^^12 is, 



{mi,m2\Ul2{Ji+ <Si l2)U^2\ni,n2) 
. e 



= exp 
= exp 
= exp 

The above expression can also be written as, 



i-(mi — ni)m2 
.J 

i-(mi — ni)m2 
. J 



e 

i-m2 
. J 



(mil Ji+|ni)(5 

^ni ^mi,ni + l^m2n2 • 



(mi,m2|C/i2(>^i+ <8i /2)t^i2|ni,n2) = (wi, m2| Ji+ (gi exp ( i ^ J^^ ) |ni, 712) 



Ul2{Ji+ O ^2)t^i2 = "^1+ exp ( z- 



Therefore, 



uI{Ji+(3I2)Ut = (C/l«)f/2)^ 

= {uIj,+u,) 



Ji+ (S) exp { i-J 

J 



[Ui (3 U2) 



C^2 exp (^^jJz2^ U2 



Now, 



U{ J{'+U(, where J^^ = U^' Ji+U^. 



fct 



In {|mi)} basis, Ji'_^_ can be written as. 



(mi|Jf+|ni) = {mi\uf Ji+U^\ni), 



= exp 
= exp 
= exp 



■ ( 2 2\ 

I— (mi -nj 



2i 

.k 



{ml - nl) 



ni 



(mi|Ji+exp 



1-^-1 + 2 



J 



(mil Ji+|ni), 

1^ 



J"_|_ = Ji_)_ exp 



i- J: 



1 



Therefore, 



t/^Ji+[/i = ufji+exp 



(A4) 



(A5) 



(A6) 



(A7) 



(A8) 



(A9) 



The operator u( is the rotation operator about y— axis with angle 7r/2, therefore U(\jxi, Jy^, Jzi)U( 
{Jzi, Jyi,—Jxi)- Hence we have 



UlJi+Ui = (J^i +iJyJexp 



J X 



(AlO) 
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Now we will calculate the other term of Ea. llA6ll . i.e., 



exp (^^^Jz2^ U2 = uCuJf exp {^'-Jz^ U^^i 



Uf exp [i^Jz.j Ui [ since [U^ J,,] = 0] 



exp -i-Jx2 
J 



since U2 is rotation matrix 



Substituting all the above results in Eq. (|A6p . we get, 

f/j,(Ji+ (g) hWr = (Jzi + iJyi)exp ij (^-Jx^ + ^ 
By taking Hermitian conjugate of the above expression, we determine 



) exp -i-Jx2 
3 



C/|(Ji_ ig) I2)Ut = exp 



i- -Jx 



1 



3 \ - 2 

Substituting, last two expressions in Ea. (jA2|) . we will get Eq.®. 



(J^i - iJy^) (g) exp ( i-Jx2 



(All) 



(A12) 



(A13) 



APPENDIX B: CALCULATION OF THE INTEGRAL PRESENT IN EQ. AND EO. (|43|l 

We know (m|z) = (to|0, 0) and using Ea. p2|l . the above mentioned integral becomes 
2j + l 



47r 



2j 

i - i 



2j \ 2j \ 2j 
j ~k J \ j -I J \ j -m 



J0=-7r 



1 + tan^ ■ 



-4i 



tan ■ 



4j — i—k — l — m 



exp[—i(j){{i + I) — (k + rn)}] sin 
After performing the 0-integral, we get 



(Bl) 



2j 



(2j + 1) 
Substituting 77 = 0/2, we get 

2(2i + 1 



2.7 \ / 2j 

J-k [j-l 



2j 

j — m 



4j+2(i+0 + l 



<5i+;,fe+m / ( cos - ) ( sin - ) dO. (B2) 



4j-2(i+i) + l 




2j \ ( 2j 



.2^' )<5.+,.+,„/"(sinr;)^^-2(^+')+i 
' / J 77=0 



(cosr;)''^+2(^+')+id?7. (B3) 



The above integral is a /3-integral, and therefore we get. 
(2j + l) 




From the relation, n) — ^(m')V{n)\ /T{m + n), we get 
2j + l 



/3[{(2j + 1) - + /)}, {(2j + 1) + + 0}]'5.+i,fc+. 



(B4) 



2j 



2j 



r(4j + 2) y V i - j / \j-k J \j -I J \j -m 



2j 



2j 



r{(2j + 1) - + 0}r{(2j + 1) + + /)}<5.+/. 



(B5) 



We know that, for any integer to, r(TO + 1) = ml. Using this relation the above expression will be equal to Ea. H35|l 
and Ea. lfl^ . 
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APPENDIX C: CALCULATION OF EQ.ll^ 



Let us first calculate the sum present in the expression of p{e). That sum can be simplified in the following way. 



exp I — i-TOir7i2 

m\.m2 — —j 

mi — 1 7712 — 1 



7V2 ^2 



exp —i-mim2 

i 



E E 



mi=— J m2 = l 



exp —1-17117712 
j 



E E 

mi — 1 m2 — —j 



exp — i-TOir7i2 

j 



E E 



-i—mim2 
j 



2N -1 



E expfi- 



2iV- 1 



.mi ,m2 — 1 



:mim2 



E 



mi ,m2 — 1 



exp —t-mim2 
j 



iV2 Ar2 



—Re ^ exp i-mi 



mi ,m2 — 1 



1^2 



(CI) 



where 'Re' denoting the real part. Now we define x = mi/ j, y = m2/ j and 5 = 1/ j, where 5 — > in large j limit. We 
can convert the above sum into an integral in the large j-limit as, 



7 lim 



x—S J y—S 

^ sin(jea;) 



dxdy cos(jexy) 
dx 



;Si(je) 



In the large j-limit, = 2j + 1 ~ 2j, therefore 



f \ 2jV-l , 2 



Ne 



If we neglect N ^ term, then we get, 



N 



Si(^ 
1 + — ^ 



(C2) 



(C3) 



(C4) 



Let us now calculate the bigger sum [see Eq. (|61ll ]. If we define li = mi — rii and I2 = m2 — n,2, then this sum will 
become, 

+M . s 

iN-\li\)iN-\l2\)exp(~ijlil2j ; M = 2j = N - 1 



hA2 = -AI 

+ M 



-M 



+M 



2N J2 iN-M+ E E {N-MiN-Mexpi-i-hh 



h = -M 



h = -M I2 = -M 



M 



4iV2M + 4Re ^ {N ~ li){N - h) exp (i-hl 



M 



h,h = l 



m'^M + 4:N^Re expli-hh) -8NRe ^ hexpli-kh 



M 



M 



4Re Z1Z2 exp ( 1-^1^2 



/l,i2 = l 



(C5) 
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We can write the first sum of the above expression as, 

^ exp ^-^1/2 j = X! ''^P ('l^^i'O ■ ^^^^ 

This sum is similar to the sum which we have calculated to derive p(e), see Eq. (|Cip . Therefore, using this previous 
result, we get the above sum as, 

^ exp (i^hh) ^ — Si(2Me). (C7) 



Now 



ll,l2 = l 



M 

e , 



Second sum = Re ^lexpl 1-/1^2 



f-i 1-1 



J 



lim / / X cos{2 Mexy)dxdy 

^-^"Jx^S Jy=S 

/ sm{2Al€x)dx 

2e Jo 



^[l-cos(2Afe)] (C8) 



M 



Third sum = Re Z1/2 exp ( 1-^1/2 



ll,l2 = l 

f-1 pi 



J 



Af ^ lim / / dxdyxy cos{2 Alexy) 

Jx=S Jy=S 

■ r N , r cos(2Mex) - 1 

/ sin 2Mea;)dx+ / , ^ dx 

Jx=s ^ ' Jx=s "^Mex 



lim 

2e <5-»o 

4^ 



— [1 - cos(2Afe) + Ci(2A/e) - ln(2Me) - 7] . (C9) 



For large j-hmit, M ~ N and therefore substituting above results in Ea. (|C5|l . we will arrive at Ea. l|?)2|l . 
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